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Abstract 



We present THC: a new high-order flux- vector- splitting code for Newtonian and special-relativistic hydrodynamics designed for direct 
numerical simulations of turbulent flows. Our code implements a variety of different reconstruction algorithms, such as the popular 
weighted essentially non oscillatory and monotonicity-preserving schemes, or the more specialised bandwidth-optimised WENO 
scheme that has been specifically designed for the study of compressible turbulence. 

We show the first systematic comparison of these schemes in Newtonian physics as well as for special-relativistic flows. In particular 
we will present the results obtained in simulations of grid-aligned and oblique shock waves and nonlinear, large-amplitude, smooth 
adiabatic waves. We will also discuss the results obtained in classical benchmarks such as the double-Mach shock reflection test in 
Newtonian physics or the linear and nonlinear development of the relativistic Kelvin-Helmholtz instability in two and three dimen- 
sions. Finally, we study the turbulent flow induced by the Kelvin-Helmholtz instability and we show that our code is able to obtain 
well-converged velocity spectra, from which we benchmark the effective resolution of the different schemes. 
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1. Introduction 

Numerical relativistic hydrodynamics has come a long way since 
the pioneering works by May & White (1966) and Wilson (1972) 
and it is now playing a central role in modelling of systems in- 
volving strong gravity and/or flows with high Lorentz factors. 
Example of applications are relativistic jets, core-collapse super- 
novae, merger of compact binaries and the study of gamma-ray 
bursts [see Marti & Muller (2003) and Font (2008) for a com- 
plete overview] . 

In all of these areas the progress has been continuous over 
the past few years to the point that relativistic computational 
fluid dynamics is starting to provide a realistic description of 
many relativistic-astrophysics scenarios [see, e.g. Rezzolla et al. 
(2011)]. Key factors in this progress have been the switch to 
more advanced and accurate numerical schemes, and in partic- 
ular the adoption of high resolution shock capturing (HRSC) 
schemes (Marti et al. 1991; Schneider et al. 1993; Banyuls et al. 
1997; Donat et al. 1998; Aloy et al. 1999; Baiotti et al. 2003) 
and the progressive inclusion of more "physics" for a more ac- 
curate description of the different scenarios. Examples of the lat- 
ter are the inclusion of magnetic fields (Koide et al. 1999; Del 
Zanna et al. 2003; Gammie et al. 2003; Komissarov 2004; Duez 
et al. 2005; Neilsen et al. 2006; Anton et al. 2006; Giacomazzo 
& Rezzolla 2007) the use of realistic tabulated equations of state, 
[see, e.g. Sekiguchi et al. (2011)], and the description of radia- 
tive processes (Farris et al. 2008; Sekiguchi 2010; Zanotti et al. 
2011). 

We expect that both improved physical models and better 
numerical techniques will be key elements in the future gener- 
ation of codes for relativistic astrophysics. On the one hand it 
is necessary to take into account many physical phenomena that 
are currently oversimplified and, on the other hand, higher ac- 



curacy is necessary to make quantitative predictions even in the 
case where simplified models are used to describe the objects 
of study. For example, in the case of inspiralling binary neutron 
stars, waveforms that are sufficiently accurate for gravitational- 
waves templates are just now becoming available and only in 
the simple case of poly tropic stars (Baiotti et al. 2009, 2010; 
Bernuzzi et al. 2012b). Clearly, even higher accuracy will be 
required as more realistic equations of state are considered or 
better characterisations of the tidal effects are explored (Baiotti 
et al. 201 1 ; Bernuzzi et al. 2012a). 

For this reason the development of more accurate nu- 
merical tools for relativistic hydrodynamics is an active and 
lively field of research. Most of the effort has been di- 
rected towards the development of high-order finite-volume 
(Tchekhovskoy et al. 2007; Duff ell & MacFadyen 2011) and 
finite-difference (Zhang & MacFadyen 2006; Del Zanna et al. 
2007) schemes, but many alternative approaches have been 
also proposed, including finite-element methods (Mann 1985; 
Meier 1999), spectral methods (Gourgoulhon 1991), smoothed- 
particle-hydrodynamics (Siegler & Riffert 2000; Rosswog 2010) 
and discontinuous Galerkin methods (Dumbser & Zanotti 2009; 
Radice & Rezzolla 2011). 

The use of flux-conservative finite-difference HRSC 
schemes is probably the easiest way of increasing the (for- 
mal) order of accuracy of the current generation of numerical 
codes: finite-difference schemes are much cheaper then high- 
order finite- volume codes since they do not require the solution 
of multiple Riemann problems at the interface between differ- 
ent regions (Shu 1999, 2001) and they are free from the com- 
plicated averaging and de-averaging procedures of high-order 
finite-volume codes [see, e.g. Tchekhovskoy et al. (2007)]. 

Here we present a new code, the Templated-Hydrodynamics 
Code (THC), developed using the Cactus framework (Goodale 
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et al. 2003), that follows this approach. THC employs a state-of- 
the-art flux- vector splitting scheme: it uses up to seventh-order 
reconstruction in characteristics fields and the Roe flux split 
with a novel entropy-fix prescription. The "templated" aspect 
reflects the fact that the code design is based on a modern C++ 
paradigm called template metaprogramming, in which part of 
the code is generated at compile time. Using this particular pro- 
gramming technique it is possible to construct object-oriented, 
highly modular, codes without the extra computational costs as- 
sociated with classical polymorphism, because, in the templated 
case, polymorphism is resolved at compile time allowing the 
compiler to inline all the relevant function calls, see e.g. Yang 
(2000). Among the different reconstruction schemes that we im- 
plemented are the classical monotonicity-preserving (MP) MP5 
scheme (Suresh & Huynh 1997; Mignone et al. 2010), the 
weighted essentially non oscillatory (WENO) schemes WEN05 
and WEN07 (Liu et al. 1994; Jiang & Shu 1996; Shu 1997) 
and two bandwidth-optimized WENO schemes: WEN03B and 
WEN04B (Martin et al. 2006; Taylor et al. 2007), designed for 
direct simulations of compressible turbulence (we recall that the 
number associated to the different methods indicates the putative 
order of accuracy). 

By far the largest motivation behind the development of THC 
is the study of the statistical properties of relativistic turbulence 
and the determination of possible new and non-classical features. 
In a recent paper (Radice & Rezzolla, in prep.), we have pre- 
sented the results of direct numerical simulations of driven tur- 
bulence in an ultrarelativistic hot plasma obtained using THC. 
More specifically, we have studied the statistical properties of 
flows with average Mach number ranging from ~ 0.4 to ~ 1.7 
and with average Lorentz factors up to ~ 1.7, finding that flow 
quantities, such as the energy density or the local Lorentz fac- 
tor, show large spatial variance even in the subsonic case as 
compressibility is enhanced by relativistic effects. We have also 
uncovered that the velocity field is highly intermittent, but its 
power- spectrum is found to be in good agreement with the pre- 
dictions of the classical theory of Kolmogorov. Overall, the re- 
sults presented in Radice & Rezzolla in prep, indicate that rela- 
tivistic effects are able to enhance significantly the intermittency 
of the flow and affect the high-order statistics of the velocity 
field, while leaving unchanged the low-order statistics, which 
appear to be universal and in good agreement with the classical 
Kolmogorov theory. 

In this paper we give the details of the algorithms used in 
THC and employed in Radice & Rezzolla in prep., presenting a 
systematic comparison between the results obtained using the 
above mentioned reconstruction schemes, with emphasis for the 
application of these schemes for direct simulations of relativistic 
turbulence. To our knowledge this is the first time that such a 
comparison has been done in the relativistic case. 

The rest of this paper is organised as follows. In Section 2 we 
present THC code in more detail: we discuss the numerical algo- 
rithms it uses and recall the equations of Newtonian and special- 
relativistic hydrodynamics. The results obtained with our code 
in a representative number of test cases for the Newtonian and 
special relativistic hydrodynamics are presented in Section 3. In 
Section 4 we present the application of our code to the study of 
the linear and nonlinear development of the relativistic Kelvin- 
Helmholtz instability (KHI) in three dimensions as a nontriv- 
ial application for our code and a stringent test of its accuracy. 
Finally Section 5 is dedicated to the summary and the conclu- 
sions. 



2. The Templated-Hydrodynamics Code 

We here briefly outline the numerical infrastructure adopted by 
our templated-hydrodynamics code and report the formulations 
of the equations of Newtonian and special-relativistic hydrody- 
namics we actually solve. 



2.1. Finite-differencing high-resolution shock-capturing 
schemes 

We can only provide here a minimal description of the numer- 
ical schemes used by THC and we refer the interested reader to 
Shu (1999) for a more detailed description of finite-difference 
ENO/WENO HRSC schemes, and to Mignone et al. (2010) for 
a detailed description of the finite-difference MP5 scheme. We 
also note that an approach very similar to ours was already 
adopted by Zhang & MacFadyen (2006) in the construction of 
the RAM code. 

Both in the case of the Newtonian and in that of the 
relativistic-hydrodynamic equations, we consider a system of 
hyperbolic balance-laws in the form 



dF°(u) 

at 



1=1 



dF\u) 

dx l 



(1) 



where the components of u are the so-called primitive variables, 
a set of N independent physical quantities, N being the number 
of equations composing (1), while F°(u) is referred to as the 
vector of conserved variables (or state vector), even though they 
are not strictly conserved in the presence of a source term S(u). 

As customary with grid-based codes, we solve numerically 
the set of equations (1) on a computational grid defined in a 
Cartesian coordinate system at uniformly distributed spatial po- 
sitions 



Xijjc = 0'A 1 , j'A 2 , kA 3 ) , i, j,k£Z, 



(2) 



and rewrite it using the method of lines, in a semi-discrete, di- 
mensionally unsplit form as 

dF lj,k F \-\l2,j,k~ F )+\l2,j,k 
~ bi,j,k + 



dt 



A 1 

r i,j-\/2,k r i,j+\/2,k r i,j,k-\/2 r ij^+1/2 



A 2 



A 3 



(3) 



where if/^ is the value of a generic quantity, if/, at while 
(F]_ l/2jk - F] +l/2jk )/A l is an high-order, non-oscillatory 1 , ap- 
proximation of -dF l /dx l at Xij^, to be specified. We use 
the method of lines to evolve (3) using a third-order strong 
stability-preserving Runge-Kutta (SSP-RK) scheme (Gottlieb 
etal. 2009). 

To illustrate how we compute the discrete derivatives in the 
right-hand- side of (3) it is useful to make a step back and con- 
sider, first a simpler scalar hyperbolic equation in one dimension, 
i.e. 

du df(u) 



dx 



= 0. 



(4) 



1 We recall that a numerical method is said to be non-oscillatory if 
N(u n+1 ) < N(u n ), where N(u n ) denotes the number of local extrema 
of some discrete representation of the function u at time t = nAt. An 
essentially non-oscillatory method is a method that is non-oscillatory on 
average only. Note that here the non-oscillatory feature is really referred 
to the reconstructed derivative of u. 
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Introducing a uniform grid, x t = /A, and define, for any function, 
v(x) the volume averages 



I rXi+i/2 
Vf = — 



v(x) dx. 



(5) 



A reconstruction operator, 7?, is a nonlinear operator yield- 
ing an high-order approximation of v at a given point x using its 
volume averages, v t . Since v(x) can be discontinuous, we distin- 
guish two different reconstruction operators, K~ and called 
the left-biased and right-biased reconstruction operators, such 
that 



pT({v/})](*)= lim v(y) + 0(A r ), 
[^ + ({v;})]«= lim v(y) + 0(A r ), 

y^>x + 



(6) 
(7) 



where we have indicated with ({vj) the notion that <R. is an 
operator that acts on a set of averages v i9 and where r is the or- 
der of the reconstruction operator <R. Hereafter we will use the 
notation vT + 1/2 and v+_ 1/2 to denote the reconstructed values in 
Xi+i/2 using K~ and ?? + respectively. Example of such operators 
are given by the ENO/WENO and the MP5 algorithms. 

In classical WENO schemes, the reconstruction is obtained 
from a weighted average of a set of lower-order reconstructions 
of the v;'s on a number of overlapping stencils. The weights are 
computed using some nonlinear smoothness indicators designed 
in such a way that the maximum possible order of accuracy is 
obtained in the case of smooth solutions. On the other hand, 
when discontinuities are detected, the order is automatically re- 
duced to avoid spurious oscillations. The bandwidth-optimized 
WENO schemes differ from the classical ones because they use 
a more symmetric stencil and because their weights are not de- 
signed to yield the maximum possible formal order of accuracy 
for smooth solutions, but are instead tuned to minimise the at- 
tenuation of high-frequency modes. In other words, while in the 
classical WENO case the weights for a smooth function are cho- 
sen to match as many terms as possible in the Taylor expan- 
sion of the target function, in the bandwidth-optimized case the 
coefficients are chosen to yield the best possible approximation 
of the Fourier coefficients of the function to be reconstructed. 
The nonlinear smoothness indicators are also modified to avoid 
"over-adaptation" of the scheme and minimise the amount of nu- 
merical dissipation [see Martin et al. (2006); Taylor et al. (2007) 
for details]. 

The MP5 scheme, on the other hand, is based on a fifth-order 
reconstruction combined with a flattening procedure designed 
to avoid the creation of artificial extrema in the function to be 
reconstructed. The monotonicity-preserving reconstruction has 
the nice property that no spurious oscillations can be produced 
by the reconstruction, which is not guaranteed for the WENO 
schemes. On the other hand the limiting procedure employed by 
MP5 requires the use of a series of conditional statements in the 
code, that are not present in the WENO case. 

The reconstruction operators are the core ingredients of both 
finite- volume and finite-difference schemes. In a finite- volume 
scheme they are used to compute the left and right state to be 
used in the (usually approximate) Riemann solver to compute 
the fluxes. In a finite-difference scheme, instead, they are used to 
compute the above mentioned non-oscillatory approximation of 
df/dx. Following Shu & Osher (1988) we introduce a function 
h(x) and such that 



/[■ 



I rXi+i/2 

A JXi-i/2 



(8) 



that is, the average of h(x) between x,-_i/2 and x;+i/2 corresponds 
to the value of / at x ( . Next, we note that 



df 
dx 



h(x i+ i /2 ) - h(xi-i /2 ) 



(9) 



where both (8) and (9) are exact expressions. Hence, by using 
the usual reconstruction operators H of order r to reconstruct 
one obtains a corresponding accurate approximation of 
order r of the derivative df/dx at jc,-. Note that h is never actually 
computed at any time during the calculation as we only need the 
values of / at the gridpoints, i.e. f[u(xi)]. 

In order to ensure the stability of the resulting scheme, one 
has to take care to upwind the reconstruction appropriately. Let 
us first consider the case in which f'(u) = df/du > 0. If we set 



'Xi+l/2 



l r 

Vi = f[u(Xi)] = - 

A Jxi-i/: 

and 

fi+l/2 = Vj +l/2 , 

then 

df(u) _ fi+l/2 ~ fi-l/2 

dx A 



fi-l/2 = V;_ 1/2 , 



+ 0(A r ), 



(10) 



(11) 



(12) 



gives the wanted high-order approximation of df/dx at X[. 

In the more general case, where the sign of f\u) is undeter- 
mined, in order to compute f+1/2, we have to split / in a right- 
going, / + , and a left-going, /", flux, / = / + +/", and use the ap- 
propriate upwind-biased reconstruction operators separately on 
both parts, in order to guarantee the stability of the method. 

There are several different ways to perform such a split and 
in our code we implemented two of them: the Roe flux- split, i.e. 



f = f\ if[/m, 1/2 ^o, 



(13) 



where Wj+1/2 = \{u t + u i+ \), and the Lax-Friedrichs or Rusanov 
flux-split (Shu 1997), i.e. 



j* = f(u) ± au, a = max[/'(w)] . 



(14) 



where the maximum is taken over the stencil of the reconstruc- 
tion operator. The Roe flux-split is less dissipative and yields a 
computationally less-expensive scheme, since only one recon- 
struction is required instead of two, but its use can result in the 
creation of entropy-violating shocks in the presence of transonic 
rarefaction waves [see, e.g. Leveque (1992)], and it is also sus- 
ceptible to the carbuncle (or odd-even decoupling) phenomenon 
(Quirk 1994). To avoid these drawbacks, we switch from the Roe 
to the Lax-Friedrichs flux split when u or / are not monotonic 
within the reconstruction stencil. 

We note that the condition that we use to switch from the Roe 
to the Lax-Friedrichs flux split is weaker than the commonly em- 
ployed condition on the sign of f'(u), [see, e.g. Leveque (1992)], 
in the sense that it results in a more frequent use of the Lax- 
Friedrichs split with respect to the usual one. According to our 
experience this prescription works very well: it is computation- 
ally less expensive to compute with respect to the standard one, 
since u and / are already evaluated on the grid, while f'(u) is 
not, and it seems to be sufficient to avoid the carbuncle phe- 
nomenon in all the tests that we performed. All the results that 
we are going to present in this paper have been obtained using 
this Roe- split with this "entropy fix". 
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We now go back to the more general system of equations 
(1). The derivatives dF c }j k /dx a can be computed following the 
procedure outlined above on a component-by-component basis. 
This approach is commonly adopted in the case of low-order 
schemes, but it often results in spurious numerical oscillations in 
the high-order (usually higher then second) case. To avoid this 
issue the reconstruction should be performed on the local char- 
acteristic variables of the systems. To avoid an excessively com- 
plex notation, let us concentrate on the fluxes in the x-direction; 
in this case, to reconstruct F\ + 1/2 . k we introduce the Jacobian 
matrices 



A a = 



dF a 
du 



a = 0,1, 



where 



(15) 



(16) 



is the average state at the point where the reconstruction is to 
be performed. We point out that the average (16), and appearing 
in (13) and (15), is much simpler than the average of u^k and 
Ui+ijfi suggested by Roe (1981). Zhang & MacFadyen (2006) 
have checked that the use of (16) in place of the average sug- 
gested by Roe (1981) does not influence significantly the quality 
of the solution in the case of finite-differencing schemes, even in 
the relativistic case. 

Hyperbolicity of (1) implies that A is invertible and the gen- 
eralised eigenvalue problem 



[A l -A (I) A°]r (I) = 0, 



(17) 



has only real eigenvalues, A^, and N independent, real right- 
eigenvectors, #•(/) [see, e.g. Anile (1990)]. We will denote with R 
the matrix of right eigenvectors, i.e. 



R 1 -r l 



(18) 



and with L its inverse. We define the local characteristic vari- 
ables 



w = Lu , 



Q = LF\ 



(19) 



and compute Q i+ i/2j,k doing a component- wise reconstruction, 
where w is used in place of u and Q in place of / in the (14). 
Finally we set 



Fi 



i+\/2,j,k 



i+\/2,j,k • 



(20) 



This procedure is repeated in the other directions and yields the 
wanted approximations of the dF a /dx a terms in x^- All the 
results that we will be presented have been obtained performing 
the reconstruction of the local characteristic variables. 



2.2. Newtonian Hydrodynamics 

The equations of classical (i.e. Newtonian) hydrodynamics de- 
scribe the conservation of mass, momentum and energy for a 
perfect fluid. They can be written in the form (1) with primitive 
variables, 



U = \p, V, 6] , 



(21) 



where p is the density, v l the velocity and e the specific internal 
energy. The conserved variables are 



F\u) = [p, pv, E] = [p, pv, p{\ v 2 + 6)] , 



(22) 



the sources are zero and the fluxes are 
F\u) = [pv\ pVv + pS\ v\E + p)] , 



(23) 



where p is the pressure and [6 i y = S 1 ^ is the Kronecker symbol. 
The system of equations is then closed by an equation of state 
p = p(p, e) and we adopt that of an ideal fluid (or Gamma law) 



p = (T-l)pe, 



(24) 



where T is the adiabatic index of the fluid. 

The Jacobians and their spectral decomposition for the equa- 
tions of Newtonian hydrodynamics and for a generic equation of 
state, can be found, for example, in Kulikovskii et al. (2001). 

2.3. Special-relativistic hydrodynamics 

In the case of the relativistic-hydrodynamic equations it is con- 
venient to work using a system of units in which c = 1 and 
we adopt the standard convention for the summation over re- 
peated indices. We consider a perfect fluid having 4- velocity 
u = (W, Wv), with W = (1 - v'v0~ 1/2 being the Lorentz fac- 
tor. Then the rest-mass current 4-vector is given by 



J = pu, 



(25) 



where p is here the rest-mass density. The stress-energy tensor is 
given by 



T = phu ®u + pg , 



(26) 



where h = 1 + e + p/p is the specific enthalpy and g is the space- 
time metric, which we take to be that of flat spacetime, i.e. where 
the only nonzero components are the diagonal ones and given by 
gfiv = (-1,1,1,1). Because our main interest with THC is of de- 
termining the statistical properties of special-relativistic turbu- 
lence and of unveiling novel and non-classical features (Radice 
& Rezzolla in prep.) we will consider the fluid not to affect the 
spacetime geometry, which we will consider to be that of a flat 
spacetime at all times. 

Conservation of rest-mass, momentum and energy are ex- 
pressed by the vanishing of the 4-divergence of / and T 



V/ = 0, v-r = o, 



(27) 



where V is the covariant derivative associated with g. Also in 
this case, the relativistic-hydrodynamic equations (27) can be 
cast in the form (1) with primitive variables 



U = [p, V, 6] . 

The conservative variables are 

F°(u) = [D, s, t] , 

where 



(28) 



(29) 



D = pW, s=phW 2 v, r = phW 2 -p-pW, (30) 



and the fluxes are given by 

F\u) = [Dv\ sV + pS\ s l - ZV] . 



(31) 



Finally, because we are working in special relativity, sources 
in the system (27) are zero and the equations are closed by an 
equation of state, which we take again to be the ideal-fluid equa- 
tion of state (24). 

An important difference between the Newtonian and the 
(special-)relativistic hydrodynamic equations is that in the latter 
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Figure 1. Density (left panels) velocity (middle panels) and pressure (right panels) for the Newtonian strong shock test. We show both the analytic 
solution (solid line) and the numerical solutions obtained with different numerical schemes (dots). The resolution is A 1 = 1/100 and the timestep 
is A = 0.002 for all the runs. Note that the numerical solution is not down-sampled and corresponds to a rather coarse resolution. 



case there is no simple analytic expression for the inverse trans- 
formation F° — > «, leading to the primitive variables from the 
conserved ones. For this reason, we use a numerical root-finding 
procedure to recover the primitive variables from the conser- 
vatives. In particular we follow the strategy of Kastaun (2006, 
2007). The primitive variables can be easily written as function 
of the conservative variables once a value for the enthalpy, h, 
is assumed. At the same time the enthalpy can be expressed as a 
function of the rest-mass density and of the internal energy using 
the equation of state. Thus we can construct the function 



3. Numerical tests 

This Section is dedicated to the presentation of some of the re- 
sults obtained with THC in a series of tests in Newtonian and 
special-relativistic hydrodynamics. 



3.1. Newtonian hydrodynamics 

We begin with a series of tests in classical hydrodynamics, be- 
fore switching to the special-relativistic case. 



g(h) = h EO s(p(h),m)-h, (32) 

and use a one-dimensional root-finding procedure to find a self- 
consistent value of the enthalpy, since g(h) = if and only if h is 
the physical enthalpy. In particular, the root-finding algorithm 
uses a combination of the Newton-Raphson method with the 
regula-falsi and the bisection schemes: we check the Newton- 
Raphson method for convergence and, in case of failure, we 
switch to the regula-falsi. The bisection scheme is used as a 
"fail-safe" root-finder in situations where the regula-falsi is con- 
verging too slowly: this is necessary only when the values of 
the conservative variables are close to an unphysical region. In 
the large majority of cases, the Newton-Raphson method usu- 
ally converges to the required level of accuracy with only three 
iterations on average. 

The Jacobians and their spectral decomposition for the equa- 
tions of special-relativistic hydrodynamics and for a generic 
equation of state, can be found in Donat et al. (1998). 



3.1.1. Strong shock 

The first test is a classical one-dimensional shock tube: the initial 
data describes two regions filled with a T = 5/3 ideal-fluid in 
equilibrium separated by a membrane. At t = the membrane is 
removed and the two regions start to interact. The initial left and 
right states are 

(Pl, v L , p L ) = (10, 0, 10) , (p R , v R , p R ) = (1, 0, 10" 5 ) , (33) 

for t > the analytic solution consists in a right-going shock 
wave, followed by a right-going contact discontinuity and a tran- 
sonic left-going rarefaction wave. 

In Figure 1 we show with a solid line the analytic solution, 
while filled circles are used to represent the solution obtained 
at time t = 0.2 with the different numerical schemes. The grid 
resolution is A 1 = 1/100 and the timestep A = 0.002 for all 
the runs. Note that the numerical solution is not down-sampled 
and hence it corresponds to a genuinely coarse resolution. We 
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do not show the results obtained with the bandwidth-optimised 
schemes since they are basically indistinguishable from the ones 
obtained using their traditional counterpart, i.e. the solution ob- 
tained with WEN03B is basically on top of the one obtained 
with WEN05 and the solution obtained with WEN04B is iden- 
tical, at the plot scale, with the one obtained with WEN07. Since 
the results obtained with the bandwidth-optimized schemes are 
found to be very close to the ones obtained with the standard 
WENO schemes in all the shock-tube test that we performed, we 
will show only the numerical solutions obtained with WEN05, 
WEN07 and MP5 in all the cases. 

Even at this fairly low resolution, all the schemes are able 
to capture well both the shock wave and the rarefaction wave, 
showing the good behaviour of our entropy fix. The shock wave 
is captured within three gridpoints. There are no appreciable 
post-shock oscillations in the solution obtained with the MP5 re- 
construction, while small oscillations are present in the velocity 
field with WEN05 and, in particular, with WEN07. 

We note that the contact discontinuity is resolved, but not 
without oscillations. However, we should bear in mind that, al- 
though the density contrast is small across the contact disconti- 
nuity, this test is a severe one due to the high Mach number of 
the shock wave, i.e. M s ~ 360. 

3.1.2. Blast wave 

The second test is similar to the first one, but results in a much 
larger density contrast at the contact discontinuity. The initial 
data is given by 

(Pl,v l ,Pl) = (10" 3 ,0, 1), (p R ,v R9 p R ) = (10" 3 ,0, lO" 5 ), 

(34) 

and the adiabatic index is still T = 5/3. Also in this case, the 
analytic solution consists in a right-going shock wave, followed 
by a contact discontinuity and a left-going rarefaction wave. 

In Figure 2 we show the exact solution at time t = 0.015 
(solid line), as well as the numerical solution (filled circles) ob- 
tained using different numerical schemes. The grid resolution is 
A 1 = 0.01 (N = 100) and the timestep is A = 0.0001 for all the 
runs. The Mach number of the shock wave is M s ~ 190 and in 
this case all our schemes are free from major oscillations. The 
shock wave is resolved within two grid points, while the contact 
discontinuity is smeared over 5-6 grid points. 

The MP5 scheme is able to properly capture the constant 
state between the shock wave and the contact discontinuity, 
while the WENO schemes result in more "rounded" solutions. 
In particular both WEN05 and WEN07 overestimated the den- 
sity contrast. 

3.1.3. Rotated Sod test 

A genuinely three-dimensional shock-tube test in Newtonian hy- 
drodynamics is offered by the classical rotated Sod test (Sod 
1978). In this case, the adiabatic index is T = 1.4 and the right 
and left states are 

(Pl,v l ,Pl) = (1,0, 1), (pr,v r , Pr ) = (0.125,0,0.1), (35) 

the initial, ID data, is rotated by 45° about the z and y axes to 
yield a shock wave that is diagonal to the principal axes of the 
grid. The analytic solution consists in a left-going rarefaction 
wave and a right-going shock wave separated by a right-going 
contact discontinuity. 



In Figure 3 we show the analytic solution in the diagonal 
direction (solid line), as well as the numerical solutions (filled 
circles), at time t = 0.2. The spatial resolution is A' = 1/200 for 
/ = 1,2, 3 and the timestep is A = 0.000625 for all the runs. 
All the schemes are able to properly capture the main features of 
the solution: the discontinuities are captured within 1 or 2 grid- 
points and both WEN05 and MP5 are able to capture the plateau 
in the velocity. The solution obtained with the WEN07 scheme 
presents oscillations in the velocity after the shock wave and in 
the density and pressure fields downstream the contact discon- 
tinuity. Similar oscillations are also observed with WEN05 and 
MP5 when the resolution is halved. 

Overall, these tests demonstrate the accuracy of the di- 
mensionally unsplit approach that we use to treat the multi- 
dimensional case. 

3.1 .4. Double Mach reflection test 

A final test in Newtonian hydrodynamics is the double-Mach 
reflection test proposed by Woodward & Collela (1984). The 
initial data describes a right-going Mach 10 shock wave mak- 
ing a 60° angle with the computational grid and intersecting the 
x-axis at x = 1/6. A perfectly reflecting wall is placed along 
the x-axis in the x > 1/6 region, while the values along the 
other regions of the boundary are set to the pre and post-shocked 
values on the left/right side of the shock wave [see Woodward 
& Collela (1984) for details on the boundary conditions and the 
initial data] . 

We considered the numerical solutions obtained with 
WEN05, WEN07, MP5, WEN03B and WEN04B in the com- 
putational domain < x < 4, < y < 1. For each of these 
schemes we performed five runs with resolutions 120 x 30, 
240 x 60, 480 x 120, 960 x 240 and 1920 x 480. In Figure 4 
we show isocontours (with an equal spacing of 0.15 between 
and 3) of the rest-mass density obtained at the highest resolution 
at time t = 0.2 by the different numerical schemes. In order to 
ease the comparison with the results reported in Woodward & 
Collela (1984), we show only the region x < 3. 

Our code performs well in this test across all the reconstruc- 
tion schemes that we tried. All the discontinuities, including the 
contact discontinuity ahead of the jet region, are well resolved 
within few grid regions, even at the lowest resolution. As we in- 
crease the number of gridpoints we do not observe any sign of 
the carbuncle phenomenon and this seems to be an indication 
that our algorithm is able to introduce enough numerical dissi- 
pation to avoid the odd-even decoupling. 

At high-enough resolution it is possible to observe the devel- 
opment of instabilities upstream from the reflected shock wave, 
generating small scale structures along the contact discontinuity. 
The ability of the different codes to resolve these structures can 
be used as an indication of their numerical viscosity. In particular 
we can see how the bandwidth-optimized schemes gain with re- 
spect to their "standard" counterparts. For instance, WEN03B, 
which uses the same stencil as WEN05, yields a solution which 
is in qualitative agreement with the one obtained using WEN05 
at twice the spatial resolution. The same is also true if we com- 
pare WEN04B and WEN07. 

All things considered, we find that the best performance is 
given by the MP5 scheme. This algorithm has a computational 
cost which is comparable with the one of WEN05, as it uses 
the same stencil. Nevertheless the solution obtained with MP5 
is very close to the one obtained with the optimized WEN04B 
scheme, which, in 3D, is almost twice as expensive as WEN05. 
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Figure 2. The same as in Figure 1, but for the Newtonian blast- wave test. The resolution is A 1 = 1/100 and the timestep is A = 0.0001 for all the 
runs. 
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Figure 3. The same as in Figure 1, but for the Newtonian rotated Sod test. The resolution is A z = 1 /200 and the timestep is A = 0.000625 for all 
the runs. 
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Figure 4. Isocontours of the rest-mass density at time t = 0.2 for the Newtonian double Mach reflection test, obtained with different schemes. We 
show 20 contours levels equally spaced between and 3. The resolution is A z = 1/480 and the timestep is A = 1 /40000 for all the runs. 



3.2. Special-relativistic hydrodynamics 

In this Section we will present the results obtained in a series of 
tests in special-relativistic hydrodynamics. 



3.2.1 . Adiabatic smooth flow 

The first test we present is designed to show the accuracy of the 
code in the case of smooth solutions and hence to measure rig- 
orously the convergence order of THC for the different schemes 
implemented. This test is very similar to the one discussed by 
Zhang & MacFadyen (2006). 

We consider a one-dimensional, large-amplitude, smooth, 
wave propagating in an isentropic fluid, with polytropic equa- 
tion of state, 



P = Kp r , 



(36) 



where K = 100 and T = 5/3. The rest-mass density at t = is 
given by 



PoW = 



\ + exp [-1/(1 - x 2 /L 2 )] , if |x| < L; 



1, 



elsewhere; 



(37) 



where, differently from Zhang & MacFadyen (2006), the initial 
profile of the rest-mass density is chosen to be C°° but is not 
analytic. We have found this choice important to obtain the cor- 
rect convergence order at very high resolutions. Indeed, when 
adopting the same profile as in Zhang & MacFadyen (2006), we 
found that the jump discontinuity in the fifth derivative of the 




Figure 5. Analytic solution and numerical solution computed with 
WEN03B with 50 gridpoints for the case of an smooth wave in an 
adiabatic relativistic fluid. The figure shows both the initial data (top 
panel) and the solution at time t = 0.8 (bottom panel). The solid line 
represent the analytic solution, while the filled circles represent the nu- 
merical one. 



initial data prevents the WEN07 scheme from achieving a con- 
vergence order larger then five. Besides this small difference, 
our initial data is basically identical to the one used by Zhang 
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& MacFadyen (2006). The initial velocity is setup assuming that 
one of the two Riemann invariants (Anile 1990), 



I /1+v 
2 ln \—v 



1 



Vr-T"l Vr-i-Ce 



In 



(38) 



where c s is the sound speed, is constant in the whole region and 
so that v = if bcl > L. The other Riemann invariant, 



1 



In 



+ c s 



(39) 



is not constant, so that the initial data describes a right-going 
wave. 

The analytic solution can be easily found in Lagrangian co- 
ordinates, up to the caustic point, using the method of charac- 
teristics (Anile 1990), while its calculation in Eulerian coordi- 
nates involves the solution of a transcendental equation. For the 
purposes of our calculation, a good-enough approximation of 
the exact solution was obtained by computing it on a very fine 
Lagrangian grid (we have used 100, 000 gridpoints), and inter- 
polated on the Eulerian grid using a cubic spline interpolation. 
This solution is then used as the reference solution against which 
the numerical solutions obtained with THC have been compared. 

In our test we set L = 0.3 and use a computational grid in 
the region -0.4 < x < 2, evolving the initial data up to time 
t = 1.6, which is approximately the time when a caustic ap- 
pears and the solution becomes discontinuous. We performed 
this test using different schemes and different resolutions and 
we measured the L 1 -norm of the error against the reference so- 
lution. Differently from all the other tests, instead of the third- 
order SSP-RK scheme, we adopt here a fourth-order RK time 
integrator. The reason for this choice is that, since this test in- 
volves a smooth solution, a SSP time-integrator is not required. 
Moreover, as we will show in the following, the use of a more 
accurate time integrator enabled us to measure a convergence or- 
der of the spatial discretization which is not spoiled by the order 
in time, the only exception being the highest resolution run done 
with WEN07. The CFL factor is set to be C « 0.2. 

The initial rest-mass density profile, as well as the solution 
at time t = 0.8, which we take as a reference time for the mea- 
sure of the error, is shown in Figure 5, together with the solu- 
tion obtained with WEN03B using 50 gridpoints. As it can be 
seen from the Figure, the solution at the considered time is still 
smooth and our scheme is able to properly capture it very well 
even at this very coarse resolution. 

The L 1 -norm of the error, as measured at the reference time, 
is shown in Figure 6. The first thing one notices is that all our 
schemes approach the expected convergence order only asymp- 
totically, at very high resolution. The reason for this behaviour is 
in the "kinks" ahead and behind the pulse, where the numerical 
error is largest. These regions are "misinterpreted" as disconti- 
nuities by the shock-detection part of our schemes, unless they 
are resolved with enough gridpoints. 

The best performing scheme in this test is the MP5 one: at 
low resolution it yields a smaller error then the seventh-order 
WENO scheme, which, in turn, is able to attain an higher conver- 
gence order only at a resolution which is unfeasible in any prac- 
tical multidimensional applications. The bandwidth optimized 
schemes present a somewhat errant behaviour in their conver- 
gence order, with WEN04B, not even showing a monotone trend 
in L 1 -norm of the error as a function of the number of gridpoint. 
We do not presently have an explanation we find sufficiently con- 
vincing for the behaviour shown. 
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Figure 6. L 1 -norm of the error for different resolutions and for different 
numerical schemes for the case of a smooth simple wave in an adiabatic 
relativistic fluid 
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Figure 7. Order of convergence, as measured using the two highest res- 
olution runs, as a function of the time to the caustic for different numer- 
ical schemes for the smooth simple wave test. 



Finally it is interesting to study the convergence order of the 
different numerical schemes as a function of t - t c , where t c is 
the time when the caustic is formed. The convergence order is 
computed using the error with respect to the exact solution of 
the two highest resolution runs and is shown in Figure 7. At 
these very high resolutions, all the schemes appear to be con- 
verging at their nominal converge order away from the caustic, 
apart from WEN07 that appears to have already reached satura- 
tion. Its order of convergence increases up to almost seven close 
to the caustic, at time t - t c « -0.4, because there the error is 
dominated by the presence of a more steep gradient ahead of 
the pulse. Before that, the error from the spatial discretization is 
probably very close to the one from the time discretization (we 
recall that we use a fourth-order RK integrator), thus degrading 
the convergence order. 

Indeed, as the time of shock formation approaches, the or- 
der of the schemes decreases slowly to the first-order expected 
in the case of discontinuities. The bandwidth optimized schemes 
and, in particular, WEN04B show, again, an erratic behaviour 
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of their convergence order. The reasons for this behaviour are 
most probably the same ones behind the similar behaviour ob- 
served while looking at the error as a function of the number of 
gridpoints. 

3.2.2. Blast wave 

When contrasted with its Newtonian counterpart, one of the most 
striking features of relativistic hydrodynamics is that relativistic 
fluids can exhibit much stronger shock waves. For this reason, it 
is important to assess the capability of the code to handle very 
strong shocks. As a first example we consider a one-dimensional 
shock-tube where the initial data is given by, 

(Pl,v l ,Pl) = (10" 3 ,0, 1), (p R ,v R ,p R ) = (10" 3 ,0, lO" 5 ), 

(40) 

and the adiabatic index is still T = 5/3. This initial data is for- 
mally identical to the one used in Section 3.1.2 in the Newtonian 
case. The analytic solution consists again in a transonic left- 
going rarefaction wave and in a right-going shock wave sepa- 
rated by a right going contact discontinuity. The shock wave has 
a relativistic Mach number M s ~ 50 (Chiu 1973; Konigl 1980). 

The results obtained with the different numerical schemes 
are reported in Figure 8, where we show the analytic solution 
(solid line) as well as the numerical ones (filled circles) obtained 
with WEN05, WEN07 and MP5, at time t = 0.4. As in the 
Newtonian case, we do not show the results obtained with the 
bandwidth optimized schemes as they are basically identical to 
the ones of their traditional variant. The spatial resolution that 
we use is A 1 = 1 /400 and the CFL factor is C = 1/5 in all cases. 

The CFL factor in this test is basically constrained by the 
MP5 scheme: while the WENO schemes appear to be robust and 
stable up to CFL factor of C « 2/5, the MP5 algorithm produces 
large oscillations and yields non-physical values in the conserva- 
tive variables unless a smaller CFL factor is used. We point out 
that for this particular test the necessity of using small timesteps 
to avoid the creation of unphysical states has been observed also 
when using other numerical schemes. For instance, when using 
the finite- volume code Whisky (Baiotti et al. 2003, 2005), with 
the HLLE approximate Riemann solver [see, e.g. Toro (1999)] 
and the PPM reconstruction (Colella & Woodward 1984), the 
maximum allowed CFL factor was also found to be C ~ 2/5. 
It is also known that, for the monotonicity-preserving property 
to hold for the MP5 scheme, the timestep must be subject to 
an additional constraint which is distinct from the standard CFL 
condition (Suresh & Huynh 1997). Yet, it is somewhat surpris- 
ing to observe that MP5 requires a timestep which is smaller 
by a factor of order two with respect to the WENO schemes. 
Furthermore, this property does not seem to be a peculiarity of 
this specific problem, since we observed a similar behaviour also 
for the other tests that we performed. 

Zhang & MacFadyen (2006) report that the use of MP5 re- 
construction results in large numerical oscillations in their spe- 
cial relativistic code. In our code, we see a similar behaviour 
unless we use a timestep which is about half of the one con- 
sidered "safe" for the WENO scheme. If the timestep is suffi- 
ciently small, on the other hand, the MP5 algorithm results in 
very accurate solutions, as in the Newtonian case. In particular 
in Figure 8 we can see the results obtained with this particu- 
lar test. Our code, with MP5, is able to capture the shock wave 
within three gridpoints while the contact discontinuity is spread 
across six points. WEN07 yields a solution of similar quality, 
while WEN05 is slightly more diffusive. 



The capability of a numerical code to capture the density 
contrast for this particular test is a classical benchmark for rel- 
ativistic hydrodynamics codes. Again, using as reference the 
Whi sky code, the use of the HLLE solver with PPM reconstruc- 
tion and with artificial compression, leads to a maximum density 
which is 71% of the analytic solution. At the same resolution, 
our THC using a WEN05, WEN07 and MP5 reconstruction is 
able to attain a maximum density which is respectively 78%, 
90% and 91% of the analytical value. These results are in very 
good agreement with the ones reported by Zhang & MacFadyen 
(2006), who measured a relative value of 72% and 79% for their 
implementation of PPM and WEN05 schemes. 

3.2.3. Shock-heating 

An even more striking example of how relativistic effect can en- 
hance the density contrasts in shock waves is given by the clas- 
sical shock-heating test. In this case, the initial data is given by 

(Pl,v l , Pl ) = (10" 3 ,v,10 3 ), (p R ,v R ,p R ) = (10" 3 ,-v,10 3 ), 

(41) 

the poly tropic index is T = 4/3 and 



v = y 1 - ~ 0.99999949 , (42) 

where the Lorentz factor is set to be W = 1000. In this case, the 
analytic solution is represented by two shocks whose collision 
compresses the fluid, converting its kinetic energy into thermal 
energy, that is, through a "shock heating". 

In Newtonian hydrodynamics the maximum compression ra- 
tio can be computed as 

^newt = = 7 , (43) 

for all the values of v, while it is easy to show (Marti & Muller 
2003) that in the relativistic case the compression ratio is 

a = fry + fTT (w_ 1} " 4003 ' (44) 

thus about three orders of magnitude larger for the same adia- 
batic index and growing linearly with the Lorentz factor. 

The exact solution at time t = 0.4, as well as the numeri- 
cal solutions obtained with our code, are shown in Figure 9. As 
it can be seen from the Figure, THC is able to handle very well 
this large compression ratio. The WEN05 and WEN07 solu- 
tions are affected by some small wall-heating effect (Noh 1987; 
Rider 2000), resulting in a slight under-density near x = and 
some numerical oscillations in the pressure. The MP5 scheme, 
on the other hand, yields a solution which is essentially free from 
oscillations in the pressure and much less affected by the wall- 
heating effect in the density variable, although at the cost of a 
smaller timestep, as discussed in the previous Section. 

3.2.4. Transverse shock 

Another peculiar difficulty of relativistic hydrodynamics and 
without a Newtonian counterpart, is that the equations for the 
momentum in the different directions are coupled together by 
the Lorentz factor: even in one-dimensional problems the appli- 
cation of a transverse velocity can change completely the so- 
lution. This feature was first pointed out by Pons et al. (2000) 
and Rezzolla & Zanotti (2002), and then used by Rezzolla et al. 
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Figure 8. The same as in Figure 1, but for the relativistic blast- wave test. The resolution is A 1 = 1 /400 and the timestep is A = 0.0005 for all the 
runs. 
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Figure 9. The same as in Figure 1, but for the relativistic shock-heating test. The resolution is A 1 = 1/100 and the timestep is A = 0.001 for all 
the runs. 



11 



David Radice and Luciano Rezzolla: THC, a new high-order FD HRSC code for special relativistic hydrodynamics 



0.02 - 



WEN05 



0.01 



0.00 



WEN07 



0.00 



MP5 



0.00 



P 

I I i i | i ~~r — r 



0.02 - 



0.01 -. 



0.02 -- 



0.01 -. 



J L 



J i i i I i i i L 



J L 



I 1 1 , 1 I , 1 ' 



-0.4 0.0 



I I i i | in — r 



J i i i I i i i L 



J i i i I i i i L 



0.4 




-0.4 0.0 



0.4 0.0 



Figure 10. The same as in Figure 1, but for the relativistic blast-wave test. The resolution is A 
runs. 
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(2003) and Aloy & Rezzolla (2006) to discover new physical ef- 
fect [see also Mignone et al. (2005); Zhang & MacFadyen (2006) 
for a description of the numerical consequences of this property]. 

To explore the flow dynamics in this case, we consider the 
same initial data as for the blast- wave test [see Section 3.2.2], 
with the only difference that we add a transverse velocity to the 
initial data, i.e. 



v* =0, 



vl = 0.99 . 



(45) 



As discussed in Zhang & MacFadyen (2006) this is not a very 
challenging test, since the velocity is added only to the cold 
fluid and it does not interact with the contact discontinuity. 
Nevertheless it is a good test to evaluate the capability of the 
code to handle the presence of a transverse- velocity at a moder- 
ate resolution, while more extreme configurations require a res- 
olution which is unreasonably high for multi-dimensional appli- 
cations to obtain decent solutions. A way around this issue is the 
use of adaptive-mesh-refinement (AMR) (Zhang & MacFadyen 
2006), which however would not be useful for our use of THC 
to study of relativistic turbulence Radice & Rezzolla in prep, (it 
is in fact debatable whether the use of AMR is advantageous in 
this case). 

In Figure 10 we show the analytic and numerical solutions 
for this problem at time t = 0.4. The presence of the transverse 
velocity widens the state between the shock wave and the contact 
discontinuity. The density contrast is also smaller with respect 
to the case where no tangential velocity is present. THC is able 
to capture the solution even at the low resolution, A 1 = 1/100, 
shown in the Figure. The MP5 scheme overestimates slightly the 
density contrast, but all of the algorithms are able to capture the 
correct location of the shock wave. 



3.2.5. Spherical explosion 

As an example of a test involving non-grid-aligned shocks we 
consider the classical test of the spherical explosion in relativis- 
tic hydrodynamics. The initial data in this case is given by 



(p, v, p) = 



(1,0,1), 
(0.125,0,0.1), 



if r < 0.4 ; 
elsewhere . 



(46) 



Since no analytic solution is known in this case, we 
use as reference solution the one computed using the 
one-dimensional spherically- symmetric discontinuous-Galerkin 
code EDGES (Radice & Rezzolla 2011) using 2000 elements of 
degree 3, solid lines, and compare it with the numerical solutions 
obtained in three-dimensions with THC in the diagonal direction. 
Both solutions at time t = 0.25 are shown in Figure 11, when 
using a resolution A' = 1/100 and a timestep A = 0.00125 
for all the numerical schemes. As in the one-dimensional case, 
a small timestep is necessary in order to avoid numerical oscil- 
lations with the MP5 algorithm, while the other schemes appear 
to be stable even with a timestep which is twice as large. 

Overall, all the schemes are able to properly capture the ref- 
erence solution even at this very low resolution: the contact dis- 
continuity is captured within two grid points and no oscillations 
are found in any of the physical quantities. 



3.2.6. Kelvin-Helmholtz instability in 2D 

The last test that we present is a classical benchmark for mul- 
tidimensional hydrodynamics codes: the simulation of the KHI 
in two-dimensions, x,y. In order to ease the comparison with 
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Figure 11. Rest-mass density (left panels) velocity (middle panels) and pressure (right panels) for the relativistic spherical explosion test. We show 
the results obtained with different schemes (dotted) as well as a reference solution (solid line) obtained with the ID spherically symmetric code 
EDGES, Radice & Rezzolla (2011), using 2000 elements of degree 3. The resolution is A f = 1/100 and the timestep is A = 0.00125 for all the 
runs. 
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where a = 0.01 is the characteristic size of the shear layer and 
^shear = 0.5, corresponding to a relative Lorentz factor, i.e. the 
Lorentz factor of a part of the fluid as seen by an observer co- 
moving with the other one, of W = 2.29. The instability is seeded 
by adding a small perturbation in the transverse component of 
the velocity, 

v y (x _ f A o Khear sin(27rx) exp [- (y - 0.5) 2 /a] , if y > ; 
V X ' y ~ j-Ao Vshear sin(27rx) exp [- (y + 0.5) 2 /<r] , if y < ; 

(48) 

where A = 0.1 is the perturbation amplitude and cr = 0.1 its 
characteristic lengthscale. The adiabatic constant is T = 4/3 and 
the initial pressure is uniform, p = 1 . The rest-mass density dis- 
tribution, which is uniform in the x-direction, is instead given 
by 



Figure 12. RMS value of the y-component of the velocity during the 
linear-growth phase of the 2D relativistic Kelvin-Helmholtz instability, 
for different numerical schemes. The solid lines represent the results 
obtained at the highest resolution, 1024 x 2048, while the dashed lines 
represent the results obtained at the lowest resolution, 128 x 256. The 
timestep is chosen so that the CFL factor is C « 0.25 for all of the runs. 



the existing literature, the initial conditions are chosen follow- 
ing Beckwith & Stone (201 1). The shear velocity is given by 



X( )= ( W tanh [(y - 0.5) /a] , if y > ; 
V W \y sh ear tanh [(j + 0.5)/a] , if : 



fj<0; 



(47) 



P(y) 



Po + pi tanh [(y - 0.5) /a] , if y > ; 
po - pi tanh [(y + 0.5) /a] , if y < ; 



(49) 



where po = 0.505 and pi = 0.495, so that p = 1 in the re- 
gions with v x = 0.5 and p = 0.1 in the regions with v x = -0.5. 
Finally the computational domain is -0.5 < x < 0.5, -1 < y < 1 
and we use periodic boundary conditions in all the directions. 
Differently from Beckwith & Stone (201 1) we do not add a ran- 
dom perturbation to the initial data and we do not take into ac- 
count the effects of magnetic fields. Nevertheless, our results are 
in good agreement with the ones by Beckwith & Stone (2011) 
in the linear phase of the instability, since the magnetic field that 
they use is too weak to play a dynamical role in this phase. 



13 



David Radice and Luciano Rezzolla: THC, a new high-order FD HRSC code for special relativistic hydrodynamics 



WEN05 



WEN05-LF 



WEN03B 








0.5 -0.5 






0.5 -0.5 




Figure 13. Rest-mass density for the 2D relativistic Kelvin-Helmholtz instability test at time t = 3 and for different numerical schemes. The 
resolution is 512 x 1024 for all the schemes and the CFL factor is C « 0.25. This figure should be compared with its equivalent Figure A.l 
presented in Appendix A and at a resolution 1024 X 2048 



We performed a series of simulations using six different nu- 
merical schemes: WEN05, WEN07, MP5, the bandwidth op- 
timized WEN03B, WEN04B and using the WEN05 scheme, 
but with the Lax-Friedrichs flux- split, WEN05-LF. For each of 
these schemes, we ran with four different resolutions: 128 x 256, 
256 x 5 12, 5 12 x 1024 and 1024 x 2048. For all the runs, the CFL 
factor was taken to be C « 0.25. 

The first quantity of interest to check is the growth rate of 
the transverse velocity during the linear-growth phase of the 
KHI, as computed with the different numerical schemes. This 
can be done by comparing the evolution of the root-mean- square 
(RMS) value of the transverse component of the velocity and de- 
fined as 



(50) 



where V is the volume of the computational domain. This is 
shown in Figure 12, where, for each numerical scheme, we show 
the results taken at the lowest resolution (dashed lines) and at 
the highest one (continuous lines). First, we notice that with our 
setup the linear-growth phase of the instability lasts up to until 
t ^ 3. At that time, in fact, the transverse velocity reaches satu- 
ration and afterwards the fluid starts to become turbulent and the 
velocity shows large fluctuations. 

Mignone et al. (2009) showed the importance of including 
the contact wave in the approximate Riemann solver in the case 
of a finite- volume code. In particular, they observed that while 
the growth rate of v y was already accurate at low resolution 
when using their HLLD Riemann solver, in the cases where the 



more dissipative HLLE Riemann solver was employed the cor- 
rect growth rate was recovered only at high resolution [similar 
results were also reported by Beckwith & Stone (201 1)]. In anal- 
ogy with what is observed with finite-volume codes, we also 
remark the importance of avoiding excessive dissipation in the 
contact discontinuity by comparing the results obtained with the 
Lax-Friedrichs and the Roe flux-split when using WEN05. The 
Lax-Friedrichs flux- split underestimates the growth rate at low 
resolution and achieves a good accuracy in its measure only at 
resolution 256 x 512 (not shown in the figure). WEN05 and 
WEN07 with the Roe flux split already have a growth rate which 
is consistent with the highest resolution runs at the lowest reso- 
lution of 128 x 256. 

The behaviour of the MP5 scheme, as well as the one of the 
bandwidth-optimized WENO schemes, is more surprising: all of 
these schemes overestimate the growth of the RMS transverse 
velocity at low resolution. This problem disappears as we in- 
crease the resolution and in the 256 x 512 the growth rate is 
already consistent with the highest-resolution one for all the nu- 
merical schemes. 

Some insight about the numerical viscosity can be gained 
by looking at the topology of the flow during the linear-growth 
phase of the KHI. In particular in Figure 13 we show a colour 
map of the rest-mass density obtained with the different schemes 
at time t = 3 using the 512 x 1024 resolution. Mignone et al. 
(2009) noticed that their solution obtained with the HLLD 
Riemann solver was showing a different structure, with the de- 
velopment of secondary small-scale instabilities along the shear 
layer, and that were not present when using the more diffusive 
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HLLE solver. Similar differences were also observed in other 
works [see, e.g. Agertz et al. (2007); Springel (2010)], and the 
presence (or absence) of these secondary instabilities has been 
often used as an indication of the numerical viscosity of the 
codes. Beckwith & Stone (2011) interpret these differences as 
an indication that HLLE is converging to a different weak solu- 
tion from HLLD. Since both solvers are entropic (i.e. with non- 
decreasing entropy), this would imply the non-uniqueness of the 
entropic solution for the Euler equations in this case. 

In agreement with the conclusions of Beckwith & Stone 
(2011), we also note that these secondary instabilities, although 
only numerical artifacts (see below), appear only in schemes 
able to properly treat the initial contact discontinuity. As a re- 
sult, the rest-mass density in Figure 13 obtained with WEN05 
and WEN05-LF match very well the ones they obtained us- 
ing HLLD and HLLE, respectively. However, our results do not 
support their conjecture that different schemes are converging 
to different weak solutions of the equations. The reason is that 
these secondary instabilities appear not to be genuine features 
of the solution and, rather, tend to disappear as the resolution 
is increased. For instance, the number of secondary vortices 
seems to change in a non-predictable fashion with the differ- 
ent numerical schemes and also with the resolution. This can 
be seen in Figure 13, which shows the great variability in the 
solutions obtained with different schemes. A similar variability 
is also present in results obtained at different resolutions with 
each scheme, WEN05-LF being the only exception. Finally, we 
point out that the size of the vortices is also shrinking as the 
resolution is increased. In essence, therefore, all of these consid- 
erations lead us to the conclusion that the secondary instabili- 
ties are triggered by the nonlinear dissipation mechanism of the 
different schemes, emerge neatly when computed with numeri- 
cal schemes that treat properly the initial contact discontinuity, 
but do not have a physical meaning (see also the discussion in 
Appendix A). 

A similar interpretation is also given by McNally et al. 
(2011), who go one step further and suggest to add additional 
viscosity to numerical codes displaying these secondary instabil- 
ities in order to prevent their growth. While the addition of extra 
dissipation is going to smooth out small scales numerical pertur- 
bations, we argue that this issue can only be resolved with the 
inclusion of physical viscosity. Artificial viscosity would proba- 
bly compete with the internal, non-linear, dissipation mechanism 
of the schemes yielding results that would be even more difficult 
to interpret (as the viscous scale will now depend both on the 
resolution and on the artificial viscosity strength). 

A more quantitative way of estimating the numerical viscos- 
ity of the code in this test has been proposed by Beckwith & 
Stone (201 1) and is based on the study of one-dimensional inte- 
grated power-spectra. Given a quantity u(x, y) we define its inte- 
grated power- spectrum 



Pl(k) = f m,y)\ 2 dy, 

k is the wa 
u(k,y) = I u(x,y)e~ 

J-l/2 



(51) 



where k is the wavenumber and 

-1/2 

- 2nikx dx, (52) 

-1/2 

is the one-dimensional Fourier transform of u. To ease the com- 
parison with the spectra reported by Beckwith & Stone (2011), 
the power spectrum is then normalized so that 

N/2 

y>i(*)=i, (53) 
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Figure 14. One dimensional power spectrum of the rest-mass density 
for the 2D relativistic Kelvin-Helmholtz instability test at time t = 3 
and for different numerical schemes. The resolution is 1024 X 2048 for 
all the schemes and the CFL factor is C « 0.25. 



where N is the number of gridpoints. The one-dimensional 
power spectrum can be used to quantify the typical scale of struc- 
tures, such as the secondary vortices discussed above, stretched 
in the direction of the bulk shear flow. 

In Figure 14 we show the spectra of the rest-mass density for 
the different schemes at time t = 3 and for the highest resolution. 
As expected the WEN05-LF scheme stands out as the scheme 
having the largest dissipation. More surprising is the behaviour 
of the bandwidth optimized schemes: they appear to be improv- 
ing over their classical counterparts only at high wavenumbers, 
that is, at scales that are dominated by the (non-physical) dissi- 
pation of the scheme. Even more unexpected is the ability of the 
MP5 scheme to resolve small scales structures and that, on the 
basis of the argument about the development of the secondary in- 
stabilities, should be more dissipative than WEN04B, but which 
instead appears to yield more small-scale structures in the rest- 
mass density. 

A similar conclusion can be obtained by studying the spec- 
trum of the kinetic energy 



E k =pW(W-l), 



(54) 



but which we do not report since its behaviour is very similar to 
the one shown for the rest-mass density. 



4. The relativistic Kelvin-Helmholtz instability in 3D 

As an example of a non-trivial application of THC and a perfect 
playground to evaluate the performances of the different numer- 
ical schemes for the simulation of turbulence reported in Radice 
& Rezzolla in prep., we present here a study of the relativistic 
turbulence generated by the KHI in three-dimensions. Our anal- 
ysis is not meant to be a systematic assessment of the accuracy of 
these methods for direct numerical simulations of compressible 
relativistic turbulence, as done, for instance by Kitsionas et al. 
(2009), or Johnsen et al. (2010) in the case of classical turbu- 
lence. Rather, our analysis is meant to assess how the different 
methods reproduce the same turbulent initial-value problem and 
to provide some insight on the spectral properties of the different 
schemes. 
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Figure 15. RMS value of the v-component of the velocity during the 
linear-growth phase of the 3D relativistic Kelvin-Helmholtz instabil- 
ity, for different numerical schemes. The thick-dashed line represent 
the results obtained at the highest resolution, 1024 x 2048, in 2D with 
WEN05. The resolution is 256 x 512 x 256 for all the schemes. 



The relativistic KHI [see, e.g. Bodo et al. (2004)] is of partic- 
ular interest because of its relevance for the stability of relativis- 
tic jets [see, e.g. Perucho et al. (2007); Perucho et al. (2010)], 
and because of its potential role in the amplification of magnetic 
fields in gamma-ray bursts [see, e.g. Zhang et al. (2009)], and bi- 
nary neutron- star mergers (Baiotti et al. 2008; Giacomazzo et al. 
2009; Obergaulinger et al. 2010; Rezzolla et al. 2011). 

We consider a triply-periodic box, -0.5 < x < 0.5, -1 < y < 
1, -0.5 < z < 0.5. The initial conditions are the same ones em- 
ployed in Section 3.2.6, and the symmetry in the z-direction is 
broken with the application of random perturbations, uniformly 
distributed in the range [0,0.01], applied to the velocity in the 
z-direction v z . We simulated this system up to the time t = 30 
on a grid of 256 x 512 x 256 cells using the WEN05, WEN07, 
MP5, WEN03B and WEN04B reconstructions. In addition, we 
have also performed one run at the same resolution using the 
second-order MINMOD reconstruction. Furthermore, because 
of the high computational costs involved to check the conver- 
gence of the code, we have managed to run only one model at 
twice the reference resolution, i.e. 512 x 1024 x 512, using the 
WEN05 reconstruction scheme. 

As discussed in the previous Section, the MP5 scheme re- 
quires a timestep which is half that of the other schemes. On the 
other hand, there is no need to run the WENO schemes with such 
a small timestep and in any "real" application we will run the 
WENO scheme with the maximum possible timestep. For this 
reason, and since we are interested in the performance of the 
different schemes in their "real- world" configuration,, we have 
not used the same CFL factor for all of them, with a consider- 
able saving in computational costs. In particular, for all the runs 
we have used a CFL factor C « 0.25, with the exception of the 
one with the MP5 scheme, where we have used C « 0.125. 



4.1. The linear evolution of the instability 

First of all, we consider the evolution of the instability during its 
linear-growth phase. At this stage, the velocity perturbations in 
the direction perpendicular to the shearing one is growing expo- 
nentially and three-dimensional effects are still negligible. 



The evolution of the RMS value of the j-component of the 
velocity is reported in Figure 15, where we show the results ob- 
tained with the different schemes, as well as a reference solution 
computed in two dimensions (2D) using WEN05 on 1024x2048 
grid points. As expected, all the numerical schemes, with the ex- 
ception of MINMOD, are in very good agreement with the 2D 
solution up to the end of the linear-growth phase at time t = 3, 
when 3D effects become important and turbulence starts to play 
an important role in the dynamics. It is interesting to note that 
MINMOD, which is the most dissipative of the schemes we are 
using, is actually overestimating the growth of the KHI. This 
suggests that some care should be taken when interpreting the re- 
sults from under-resolved simulations, since it is not always safe 
to assume that the high numerical viscosity of the low-resolution 
runs tends to suppress the instability, yielding a lower-bound on 
its growth. 

The rest-mass density at time t = 3 is shown in Figure 16, 
where, in analogy with the 2D case, we find that the more dis- 
sipative schemes (i.e. MINMOD), do not show any sign of sec- 
ondary instabilities apart from the seeded ones, while the least 
dissipative ones (i.e. WEN03B, WEN04B) show the emergence 
of secondary vortices. 

At this point in time, the flow is still mostly two-dimensional, 
but it is interesting to notice the effects resulting from the small 
perturbations in the vertical velocity v z . The perturbations have 
the same statistical properties in all the different models, but their 
effects are appreciably different for the different schemes, as can 
be seen from Figure 16. Although the perturbation in v z is ran- 
dom and it does not preserve any of the symmetries of the prob- 
lem, it is still symmetric, on average, with respect to the j-axis. 
For this reason, one does not expect to find a large symmetry vi- 
olation until the time when these perturbations have had enough 
time to grow and the dynamics has started to become domi- 
nated by three-dimensional effects. Yet, the optimized schemes 
WEN03B, WEN04B appear to be much more sensitive to the 
symmetry breaking than the others. The reason for this is prob- 
ably that the bandwidth optimized- algorithms appear to trigger 
smaller-scale secondary instabilities, which, in turn, are more 
easily affected by the perturbations in the vertical velocity, since 
the perturbation does not average out at small scales. 

4.2. The nonlinear evolution of the instability 

The linear-growth phase of the KHI instability ends when the 
primary vortices become unstable to secondary instabilities and 
the flow starts the transition to turbulence. At this point, three- 
dimensional effects start to become dominant and the dynamics 
is qualitatively different from the one observed in the 2D case. 

In order to better track the evolution of the fluid in this phase, 
we monitor the concentration of a passive scalar field, 0, trans- 
ported with the fluid via the advection equation 



dt ^ r)r> 



(55) 



Equation (55) is not a conservation law, so that in can not be writ- 
ten directly 2 in the form (1), but it is nevertheless a hyperbolic 



2 Equation (55) can be written in conservation form at the price of 
introducing an additional, conserved variable, = pW(p. In our study 
this complication is not necessary, as we use the tracer only as a diag- 
nostic quantity. On the other hand, in situations where, for instance, the 
tracers are used to model the chemical composition of a fluid in a react- 
ing flow, it may be important to ensure the conservation of the different 
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Figure 16. Rest-mass density in the z = plane for the 3D relativistic Kelvin-Helmholtz instability study at time t = 3 and for different numerical 
schemes. The resolution is 256 x 512 x 256 for all the schemes. 



equation of the Hamilton- Jacobi type. For this reason, in order 
to solve numerically (55) we use a class of numerical schemes 
designed for this family of equations and that can be built using 
the non-oscillatory reconstruction of the derivative introduced in 
Section 2 [see, e.g. Shu (2007)]. In particular, the semi-discrete 
form of equation (55) becomes 



d<pi,j,k 



_ x <Pi-l/2,j,k - (Pi+\/2,j,k y <Pi,j-l/2,k ~ <Pi,j+l/2,k 
dt ~ ViJ > k ~ + r " ' ~ ' + 



A 1 

1>i,j,k-l/2 - <Pi,j,k+l/2 



v i,j,k 



A 2 



A 3 



(56) 



where (<pi-i/2j,k ~ <Pi+i/2j,k)l 'A 1 is a high-order non-oscillatory 
approximation of -dcp/dx 1 at x^k obtained using an upwind- 
biased reconstruction, i.e. 



»i-l/2,j,k 



u i-l/2,j,k> 



ifv x ^0. 



(57) 



The "fluxes" in the other directions are obviously computed in 
an analogous way. 

The time evolution of the tracer, as computed with WEN05 
at the highest-resolution is shown in Figure 17. At the initial 
time, we set the scalar field to be equal to the rest-mass den- 
sity, so that the initial configuration consists in two regions with 
opposite "colour", separated by a thin transition layer. At the end 
of the linear regime, i.e. at time t 3, the tracer profile is dis- 
torted by the presence of the primary vortices as well as of the 



species and a conservative approach may be preferred [see, e.g. Plewa 
& Mueller (1999) for an example of this approach]. 



secondary ones, but there is only a marginal mixing of the two 
"phases" of the fluid. Note that, since we do not include any ex- 
plicit dissipation term in the advection equation (55), the mixing 
of the tracer happens only due to the numerical dissipation. As 
the vortices start to become unstable, the fluid starts to concen- 
trate the scalar field in thin structures and the two regions of the 
fluid start to mix around the shearing region. As the simulation 
proceeds, the width of the region effected by the mixing grows, 
till at time t 24, when there are only small patches of the fluid 
that still have a uniform tracer (colour). At the final time, t = 30, 
there are no regions in the flow that are unaffected by the mixing 
and the initial structure is mostly destroyed, even though perfect 
mixing has not been achieved yet. 

We track the evolution of the variance of the scalar field, 
which we compute as 



Var[0] = (|0 - 



<0>| 2 >, 



(58) 



and which, for t > 5, we find to be very- well fitted by a simple 
exponential law of the type 



Var[0] = Ke~ tlT . 



(59) 



The values of the fitting constant for the WEN05 scheme at the 
highest resolution are K = 0.28 and r = 14.75. The mixing 
timescale, r, exhibits only small changes between the runs at 
different resolutions, with the exception of the results obtained 
with MP5, where the timescale is r = 17.5. Hence, the total time 
of the simulation, t = 30, is roughly equivalent to two e-folding 
times for the mixing of the passive tracer. 
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Figure 17. Evolution of the concentration of the passive tracer in the z = plane for the 3D Kelvin-Helmholtz instability obtained using WEN05 
on 512 x 1024 x 512 gridpoints. 



4.2.1. Fully-Developed turbulence 

After the linear-growth phase, the flow quickly becomes turbu- 
lent. By the time we stop the simulation at t = 30, the turbulence 
is fully developed even though the flow is still not isotropic, but 
preserving some memory of the initial configuration. Because of 
the large computational costs involved in these tests (which, we 
recall, have been performed for six different methods) we have 
not carried out the simulations for longer times, assuming that 
the properties of the non-perfectly isotropic turbulence at the fi- 
nal time is sufficiently close to the fully isotropic turbulence. We 
point out that our simulation time is more then two times longer 
than the one reported in Zhang et al. (2009), where a setup sim- 
ilar to ours was used, but in the more challenging regime of rel- 
ativistic MHD. 

By far the most interesting quantity to study is the three- 
dimensional velocity power spectrum 

P v (k)=l f \v(k)\ 2 dk, (60) 

z J\k\=k 

where v(k) is the three-dimensional Fourier transform of v(jc), 

v(Jfc) = f v(x)e- 27rikx dx, (61) 
Jv 

and V is the volume of the computational region. The integral in 
(60) is computed following Eswaran & Pope (1988) as 

w=1 2ir 2 ]m]2 > (62) 

k k-l/2<\k\<k+l/2 



io- 1 



10" 



10" 



10" 



H 1 — I I I I | 



H 1 — I I I I | 




WEN05_N512 
WENO5_N1024 



10° 



10 1 



10 2 



Figure 18. Compensated power spectrum of the velocity at time t = 
30 for the 3D Kelvin-Helmholtz instability test, using WEN05 at two 
different resolutions, 256 x 512 x 256 and 512 x 1024 x 512. 



where Nk is the number of discrete wave-numbers in the shell 
k - 1/2 < \k\ < k + 1/2. For simplicity, we study the data in 
the artificially-enlarged cubic domain [-1, l] 3 and we do that by 
duplicating original data in the x and z-directions, exploiting the 
symmetry of the problem. This procedure avoids the complica- 
tions of a computational domain which does not have the same 
extent in all directions and yields an "equivalent resolution" 
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Figure 19. Compensated power spectrum of the velocity at time t = 30 for the 3D Kelvin-Helmholtz instability test, using different numerical 
schemes. In all the cases the resolution is 256 x 512 x 256. 



of 512 3 and 1024 3 points for the low and high-resolution runs 
respectively. Clearly the first few wave-numbers of the spec- 
trum are influenced by this procedure, but all the higher wave- 
numbers are essentially unaffected. 

At the final time, the flow is subsonic (At max ^ 1) and rel- 
ativistically warm (e > 1). Under these conditions, studies of 
Newtonian (Porter et al. 2002) and relativistic (Zhang et al. 2009; 
Inoue et al. 2011; Zrake & MacFadyen 2012) transonic turbu- 
lence suggest that the velocity spectrum should be consistent 
with the Kolmogorov phenomenology (Kolmogorov 1991). In 
particular, the power spectrum should scale as 

p v (k) ~ r 5/3 , (63) 

in the so called inertial range, that is at those scales where the 
fluid motion is sufficiently independent from the large-scale dy- 
namics and from the small-scale viscosity, so as to exhibit a self- 
similar universal behaviour. 

In Figure 18 we show the compensated velocity spectrum, 
i.e. k 5/3 P v (k), at time t = 30 obtained from the data of the two 
WEN05 runs. More specifically, the low-resolution spectrum is 
shifted to larger wavenumbers by a factor two and scaled assum- 
ing a /T 5/3 scaling of the spectrum. The rationale behind this 
procedure is that we are interested in the (eventual) selfsimilar 
behaviour of the spectrum and it is therefore useful to consider 
the low-resolution run as a realisation of the same flow as the 
high-resolution one, but in a smaller volume. In other words, we 
interpret the spectrum of the low-resolution run as being com- 
puted from a subsample of the data of the high-resolution one 



[see, e.g. Bodo et al. (201 1) for a more detailed discussion of the 
issue of convergence for direct simulations of turbulent flows] . 

The plot demonstrates the statistical convergence of the sim- 
ulation, with the exception of the very low-wavenumber part 
of the spectrum, where the convergence is probably spoiled by 
the symmetry. At the same time, the plot also shows that only 
the high-resolution run seems to be able to cover a sufficiently 
wide part of the inertial range to give a clear indication of the 
Kolmogorov -5/3 scaling. 

In Figure 19 we show the compensated velocity spectra at 
time t - 30 obtained with the different numerical schemes. For 
each scheme, we highlight with a grey shaded area the part of the 
spectrum that appears to be "compatible" with the -5/3 scal- 
ing inferred from the high-resolution WEN05 run. The width 
of this region, expressed in terms of wavenumbers, is also indi- 
cated on the figure. Clearly, this measure has to be taken with a 
bit of caution, since there is a certain degree of arbitrariness in 
the identification of the "inertial range"; furthermore, the judg- 
ment is also made more difficult by the fact that all of the results 
are obtained at a resolution which probably is not high-enough 
and that a convergence study is only available for the WEN05 
case. Notwithstanding these caveats, the difference between the 
various schemes is already sufficiently clear to deduce some of 
their properties despite the limited amount of data that we were 
able to generate. 

A first conclusion to be drawn is about the importance of 
the use of high-order schemes, whch is apparent if we compare 
the different spectra with the one obtained with the MINMOD 
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scheme. This second-order algorithm, in fact, yields a spectrum 
which is completely dominated by the so-called "bottleneck- 
effect" [see, e.g. Brandenburg & Nordlund (201 1); Schmidt et al. 
(2006)], i.e. by a region where the power- spectrum shows an 
excess due to viscous effects. No clear indication of an inertial 
range appears anywhere in the spectrum, with the only "flat" 
region being the middle of the bottleneck zone. This could be 
easily mistaken as the inertial range in the absence of a proper 
convergence analysis. For this reason, and as remarked many 
times already, care should be taken while interpreting the results 
obtained with low-order schemes. 

Secondly, we observe that WEN04B has an effective res- 
olution which is about 50% larger then the one from WEN05 
and 20% larger then WEN07. Given that saving a factor 1.5 in 
resolution corresponds, in 3D to a decrease of the computation 
costs by a factor five 3 , we conclude that the use of WEN04B 
over WEN05 is well justified, since WEN04B is roughly twice 
as expensive as WEN05 in 3D. On the other hand, the spectral 
resolution of WEN04B does not appear to be better than the 
one of the MP5 scheme, which has similar computational costs 
(due to the stricter CFL limitation), but which is also expected to 
have better parallel scaling because of its more compact stencil. 
Overall, MP5 shows an "inertial-range" as large as WEN04B. 
We also note that WEN03B does not seem to yield any improve- 
ments over WEN05. 

All things considered, the main differences between the 
bandwidth-optimized schemes and their traditional counter- 
parts seem to lay in the bottleneck region: WEN03B and 
WEN04B have a much less pronounced bottleneck with re- 
spect to WEN05, WEN07 and MP5. This suggests that these 
schemes should be considered especially for under-resolved tur- 
bulent flows, where the spectrum is basically entirely domi- 
nated by the dissipation. MP5, on the other hand, can be par- 
ticularly useful for those problems where one is interested in 
well-resolved quantities, such as in direct numerical simulations 
of turbulence, since the scales affected by the numerical dissi- 
pation are more easily identified by the pronounced bottleneck. 
MP5 should also be the scheme of choice in Newtonian hydro- 
dynamics, since there its timestep constraint seems to be less 
severe. 



5. Conclusions 

We have presented THC, a new multi-dimensional, finite- 
difference, high-resolution shock-capturing code for classical 
and special-relativistic hydrodynamics. THC employs up to 
seventh-order accurate reconstruction of the fluxes in local char- 
acteristic variables and the Roe flux- vector- splitting algorithm 
with a novel entropy-fix prescription. The multi-dimensional 
case is treated in a dimensionally unsplit fashion and the time in- 
tegration is done with a third-order strongly-stability -preserving 
Runge-Kutta scheme. 

We have carried out a systematic comparison of the results 
obtained with our code using five different reconstruction oper- 
ators: the classical WEN05, WEN07, MP5, as well as two spe- 
cialised bandwidth-optimized WENO schemes: WEN03B and 
WEN04B. For all schemes, we have checked their ability to 
sharply capture grid-aligned, diagonal or spherical shock waves, 
and have carried out a rigorous assessment of their accuracy in 
the case of smooth solutions. Finally, we have contrasted the per- 

3 We note that when increasing the resolution, also the timestep is 
reduced via the CLF condition. Hence, the additional cost goes like the 
fourth power of the ratio in resolutions. 



formance of the different methods in the resolution of small scale 
structures in turbulent flows. To the best of our knowledge, this 
is the first time that such a comparison has been carried over in 
the special-relativistic case. 

Among the different tests studied, some are highly nontriv- 
ial, such those involving the linear and the nonlinear phases of 
the development of the Kelvin-Helmholtz instability for a rela- 
tivistic fluid in two and three dimensions. In particular, we have 
shown the importance of using schemes that are able to properly 
capture the initial contact discontinuity in order to obtain the 
correct growth rate of the RMS transverse velocity in the linear- 
growth phase of the instability at low resolution, confirming the 
findings by Mignone et al. (2009) and Beckwith & Stone (201 1). 

When studying Kelvin-Helmholtz instability in two dimen- 
sions, we have investigated the nature of the secondary vortices 
that appear during the initial stages of the instability when us- 
ing some of the numerical schemes considered. We have then 
clarified that these vortices are not genuine features of the so- 
lution of the equations, but rather numerical artefacts that con- 
verge away with resolution. When studying Kelvin-Helmholtz 
instability in three dimensions, we have instead investigated the 
"mixing timescale" of the instability and the subsequent turbu- 
lent flow, showing that we are able to obtain a converged mea- 
sure of the velocity power spectrum, using the WEN05 scheme. 
Our data offers a clear indication that the Kolmogorov phe- 
nomenology (Kolmogorov 1991) holds also for the turbulence 
in a subsonic relativistically-warm fluid. Using the Kolmogorov 
-5/3 scaling as a reference, we have estimated the effective iner- 
tial range of the different schemes, highlighting the importance 
of using high-order schemes to study turbulent flows. The code 
has also been used in a systematic investigation of direct numer- 
ical simulations of driven turbulence in an ultrarelativistic hot 
plasma, whose results will be presented in a companion paper 
(Radice & Rezzolla in prep.). 

Finally, THC represents the first step towards the implemen- 
tation of new and high-order methods for the accurate study of 
general-relativistic problems, such as the inspiral and merger 
of binary neutron stars (Baiotti et al. 2011) and their relation 
with gamma-ray bursts (Rezzolla et al. 2011). We are in fact 
convinced that the transition towards higher-order methods is 
now a necessary and an inevitable step for a more realistic de- 
scription of the complex phenomenology associated with these 
relativistic-astrophy sics processes . 
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Appendix A: Secondary vortices in the KHI 

In Sect. 3.2.6 we have commented that we believe the secondary 
instabilities appearing the in initial stages of the KHI are not 
genuine features of the solution, but just artefacts of the (low) 
resolution. The latter acts differently on different schemes, lead- 
ing to a non-predictable change in the number and shape of these 
secondary vortices. Convincing evidence that these are indeed 
artifacts is shown in Figure A.l, which is the same as Figure 13, 
but at the higher resolution of 1024 x 2048. When comparing 
the two figures it appears clear that these secondary instabilities 
are much smaller. Interestingly, therefore, the more dissipative 
scheme WEN05 with the Lax-Friedrichs flux-split, WEN05- 
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Figure A.l. Rest-mass density for the 2D relativistic Kelvin-Helmholtz instability test at time t = 3 and for different numerical schemes. The 
resolution is 1024 x 2048 for all the schemes and the CFL factor is C « 0.25. This figure should be compared with its equivalent Figure 13 at a 
resolution 512 x 1024 



LF, seems to converge more rapidly to the correct solution (at 
least in these initial stages) than its less diffusive counterparts. 
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